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ABSTRACT 

Context. The physics of the pulsar inner magnetosphere remains poorly constrained by obser- 
vations. Although about 2000 pulsars have been discovered to date, little is known about their 
emission mechanism. Large vacuum gaps probably exist and a non-neutral plasma made of elec- 
trons in some regions and of positrons in some other regions fills space to form an electrosphere. 
Aims. The purpose of this work is to study the stability properties of the differentially rotating 
equatorial disk in the pulsar's electrosphere for which the magnetic field is assumed to be dipolar. 
In contrast to previous studies, the magnetic field is not restricted to be uniform. 
Methods. A pseudo-spectral Galerkin method using Tchebyshev polynomials expansion is de- 
veloped to compute the spectrum of the diocotron instability in a non-neutral plasma column 
confined between two cylindrically conducting walls. Moreover, the inner wall carries a given 
charge per unit length in order to account for the presence of a charged neutron star at the centre 
of the electrosphere. 

Results. We show several eigenfunctions and eigenspectra obtained for different initial density 
profiles and electromagnetic field configurations useful for laboratory plasmas. The algorithm is 
very efficient in computing the fastest growing modes. Applications to a "cylindrical" electro- 
sphere are also shown for several differential rotation profiles. It is found that the growth rates of 
the diocotron instability are of the same order of magnitude as the rotation rate. 
Conclusions. The instability develops on a very short timescale and can account for very efficient 
particle diffusion across the magnetic field lines as already claimed in a previous work. The exact 
geometry of the confined plasma, let it be a thin disk or a cylinder, does not significantly affect 
the spectrum of the diocotron instability. 
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1. INTRODUCTION 

The detailed structure of charge distribution and electric current circulation in the closed magneto- 
sphere of a pulsar remains poorly understood. Although it is often assumed that the plasma entirely 
fills the space and corotates with the neutron star, it is on the contrary very likely that it only partly 
fills it, leaving large vacuum gaps in between plasma-filled regions. The existence of such gaps in 
aligned rotators has been very clearly established by Krause-Polstorff and Michel (|1985a|[l985bl l. 
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Since then, a number of different numerical approaches to the problem have confirmed their conclu- 
sions, including some work by Rylov d 19891 1. Shibata d 19891 1. Zachariades ( 119931 ). Neukirch ( 119931 ). 
Thielheim and Wolfsteller ( |1994| l. Spitkovsky and Arons ( |2002| l and by ourselves (Petri et al. 
I2002al) . This conclusion on the existence of vacuum gaps has been reached from a self consistent 
solution of Maxwell's equations in the case of the aligned rotator. Moreover, Smith et al. ( 120011 1 
have shown by numerical modeling that an initially filled magnetosphere like the Goldreich-Julian 
model evolves by opening up large gaps and stabilises to the partially filled and partially void so- 
lution found by Krause-Polstorff and Michel ( |1985a| l and also by Petri et al. (i2002al l. The status 
of models of the pulsar magnetospheres, or electrospheres, has been recently critically reviewed 
by Michel (I2005I I. A solution with vacuum gaps has the peculiar property that those parts of the 
magnetosphere which are separated (following magnetic field lines) from the star's surface by a 
vacuum region are not corotating and suffer differential rotation. 

This rises the question of the stability of such charged plasma flow. The differential rotation 
in the equatorial non neutral disk induces the so-called diocotron and magnetron instabilities, well 
known to plasma physicists (O'Neil [19801 Davidson Tl 9901 O'Neil and Smifli |1992] l. In the inner 
parts of the magnetosphere, far from the light cylinder, the instability reduces to its electrostatic 
form, the diocotron instability. The linear development of the diocotron instability of a thin dif- 
ferentially rotating charged disk has been studied by Petri et al (I2002b 1) and shown to proceed at 
a growth rate comparable to the star's rotation rate. The non linear development of this instabil- 
ity has been studied by Petri et al (120031 ). still in the framework of an infinitely thin disk model. 
When there is no external source of charges feeding the magnetosphere, it has been found that the 
plasma remains confined in regions close to the equilibrium state as also shown by Aly (120051 ). 
When however the disk can be fed by some charge source, Petri et al (120031 ) have shown that the 
instabihty causes a cross-field transport of these charges in the equatorial disk, evolving into a net 
out-flowing flux of charges. Spitkovsky and Arons ( 120021 ) have numerically studied the problem, 
and concluded that such charge transport tends to fill the gaps with plasma. Since a filled magneto- 
sphere is stable to the diocotron instability, it seems unlikely that the instability could indeed bring 
the system to a completely filled state, even though the disk would not remain thin, as assumed by 
Petri et al. ( |2002bl 120031 ). It should be noted that the appearance of a cross-field electric current as 
a result of the diocotron instability has been observed by Pasquini and Fajans (120021 ) in laboratory 
experiments in which charged particles were continuously injected in the plasma column trapped 
in a Malmberg-Penning configuration. 

The aim of this work is to show that the fact that the differentially rotating charged disk is not 
infinitely thin does not invalidate our former conclusion (Petri et al. l2002a l) that the growth rate of 
the instability is as fast as the star's rotation rate. The situation opposite to a thin disk is that of an 
infinitely thick disk. By this we mean a plasma column, the structure of which is invariant with the 
cylindrical coordinate z and in which the magnetic field is dependent on the cylindrical coordinate 
r and oriented parallel to the z direction. Whereas the thin disk model can conveniently describe 
perturbations the horizontal size of which is much larger than the disk's thickness, the (infinitely) 
thick model would be more appropriate to perturbations the horizontal size of which is much less 
than the disk's thickness. In the electrostatic approximation the magnetic field is not altered by the 
perturbations and remains straight. From the inertia-less approximation, the charged fluid velocity 
is given by eq. ( fTTI ) which implies that /?+vxB = and then that E B = 0. As a result the perturbed 
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electric field and flow are 2-dimensional, since the electric potential perturbation is independent of 
z. It should be kept in mind that the infinitely thick disk model with a z-aligned magnetic field is in 
fact meant to represent a disk of a finite thickness embedded in a 2D potential magnetic field. Such 
2D fields, Uke 

B^B,(r,z)e, + B,(r,z)e, (1) 

as opposed to ID fields like B - Bz('")Cz, may be current free and still have a radial dependence 
of Bz('", 0) in the equatorial plane. In the ID model, the disk, though not infinitely thin, is regarded 
as not being thick enough for the radial field component Br{r,z) to play a role. The component 
B^(r,Q) is taken as an approximation to Bj{r,z)- The domain of validity of a ID model then is 
that Br(r,z) remains small in the disk, while the horizontal size of the perturbations is less than 
the disk's thickness. The real situation is in between the infinitely thin and the infinitely thick disk 
model. Dealing with a finite thickness disk would involve the consideration of the 2D magnetic field 
structure and the numerical calculation of the solution to the Poisson's equation in the appropriate 
geometry. Obtaining this solution has been the most time consuming part of the calculations in 
Petri et al. (I2003I I and things would be even worse if a 2D geometry with a thick disk were to be 
considered. Adopting a ID unperturbed geometry model considerably alleviates such difliculties 
and seems appropriate to reach a conclusion on the range in which the growth rates of the diocotron 
instabihty are to be found when the infinitely thin disk approximation is relaxed. 

In this paper we present a numerical analysis of the linear growth of the diocotron instability 
in the ID model. We plan to compute the full non linear development of the instability in the same 
geometry, by using a 2D electrostatic PIC code, a study that would also be relevant to laboratory 
setups. 

The diocotron spectrum for cylindrical geometry and uniform external magnetic field has been 
investigated analytically by Levy (1 19651 1 and numerically by Goswami et al. (|1999b who used a 
finite difference schemes. Variational techniques for toroidal non-neutral plasmas have also been 
used, see Bhattacharyya (|2000| l. 

The paper is organized as follows. In Sec.|2] we describe the initial setup of the plasma column 
in the laboratory consisting of an axially symmetric equilibrium between two conducting walls. 
The slightly different approach used for the pulsar's electrosphere is also described. In Sec.|3] we 
recall the generalised hnear eigenvalue problem satisfied by the perturbed electric potential. Next, 
in Sec. |4] the pseudo-spectral Galerkin numerical algorithm to compute the diocotron spectrum 
is presented. Finally, some typical spectra for laboratory plasma are shown in Sec.|5] We discuss 
in more detail the application to electrospheric plasmas which is the main aim of this work. The 
conclusions and the possible generalisation are presented in Sec. |6] 

2. INITIAL SETUP 

In laboratory plasmas, the motion of particles is imposed by external devices controlling the po- 
tential and electric field. It is the most common situation. However, having in mind to applied the 
algorithm to pulsar's electrosphere, it is also interesting to study the stability properties of non- 
neutral plasmas by imposing the rotation profile instead of the density profile. Therefore, we first 
describe the case encountered in the laboratory and then discussed how to apply and extend it to 
pulsars. 
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2.1. Laboratory plasma 

We consider a one species non-neutral plasma consisting of particles with mass m and charge q 
trapped between two cyUndrically conducting walls located at r = Wi and r - 'W2 > W\. The 
plasma column itself is confined between Ri > W\ and R2 < W2. This allows us to take into account 
vacuum regions between the plasma and the conducting walls. We adopt cylindrical coordinates 
denoted by (r, (f, z) and the corresponding basis vectors {e^, e^, Cz). In order to simulate the presence 
of a charged neutron star generating a radial electric field, the inner wall at W\ carries a charge Q 
per unit length such that its electric field is simply given by Maxwell-Gauss theorem : 

£wW=— (2) 

In EQr 

In the equilibrium configuration, the particle number density is n{r) and the charge density is p{r) - 
qn{r). In contrast to earlier studies, the external magnetic field, along the z-axis, is not necessarily 
uniform 

B^B,ir)e, (3) 

Strictly speaking, this magnetic field is associated with an azimuthal current because V A B ^ 0. 
However, as claimed in the introduction, there is also a (weak) radial component Eq. ([TJ which 
should compensate for this current. Let's find the conditions for which approximation Eq. (O holds, 
assuming a power law for Bz, Eq. ( |48] l. From Eq. ([T), the (^-component of the curl reads 
dB, dB^ dB, Bz 

^ - ^ = ^ + a— (4) 

dz or oz r 

With the boundary conditions Bj{r, z = 0) = 0, we found that near the equatorial plane 

B, ^-a-B, (5) 
r 

a being of order unity, the radial component of the magnetic field is negligible whenever z r. 
Therefore, our approximation Eq. (O holds in the thin disk limit (intermediate between infinitely 
thin and infinitely thick). 

The electric field is made of two parts, the first one arising from the plasma column Ej, itself, 
and the second one from the inner conducting wall E^,, Eq. (|2]l. We assume that the electric field 
induced by the plasma vanishes at r = Wi, i.e. Ep(Wi) - 0. We therefore get 

£p(r)=— r p(r')r'dr'er (6) 

The total electric field, directed along the radial direction e,, is therefore 



E^Ep+E^^ErCr (7) 

At equilibrium, the plasma is rotating between the two walls at a speed Q.{r). The particles have 
mass /«£ and charge q. In the stationary state, the centrifugal force is balanced by the Lorentz force 
such that 

^(£'rH-rQBz)+WerQ^ =0 (8) 

Introducing the cyclotron frequency by coc = ^ and the plasma frequency by cj^ = the 
angular speed of the column satisfies the quadratic equation 

Q2 -H We Q -I- — = (9) 
nif-r 
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leading to two possible solutions for Q, namely 



2 



1± 



(10) 



If the electric field is sufliciently weak, i.e. qE^ <& mgrw^, the solutions are approximated by 
D. - -(jL>c and 

The weak field limit is equivalent to a low density plasma column because it implies Wp <s; ■ The 
solution Eq. ( fTTT ) corresponds to the usual electric drift approximation (Davidson and Tsang [T984l ). 
In the remaining of this paper, we retain this assumption. Therefore, the equation of motion for 
the charged particles is replaced by the electric drift approximation. Indeed, we consider only the 
low-density non-neutral plasma case for which the diocotron regime is valid. Neglecting electron 
inertia, the electric drift approximation applies. For 2 = and a constant density profile in the 
plasma column with Wx - 0, it corresponds to a circular motion at the diocotron frequency defined 
by, see for instance Davidson ( 119901 ) 

«o = (12) 
2 o)^ 

2.2. Electrospheric plasma 

Our purpose is to study the non-neutral plasma instabilities occurring in the pulsar inner magne- 
tosphere. The most interesting feature obtained in the work by Petri et al. ( |2002a| l by constructing 
an electrospheric model was a dififerential rotation in the positively charged equatorial disk. It can 
indeed be shown that large vacuum gaps imply a significant departure from corotation. That is why 
when applying our algorithm to pulsars, it is more appropriate to assume a given rotation profile in 
the plasma Q(r), instead of densities and potential imposed by external devices. In the electric drift 
approximation, using Maxwell-Gauss law, the charge density is recovered according to ideal MHD 
by 

E^-(£lAr)AB (13) 
leading to the charge density as follows 

(14) 



-So 



r nr V nr 



r dr ^ ' dr 

Therefore, there are only 2 independent functions, namely the couple (Sz, f^) appropriate for pulsar 
electrospheric plasmas and the couple (B^,p) useful for laboratory plasmas. 



3. LINEAR ANALYSIS 

We consider two different plasma configurations. In the first case, the plasma column is supposed 
to be in contact with both the inner and the outer wall. This will be called single domain. In the 
second case, vacuum regions exist between the plasma column and the inner and/or the outer walls. 
This will be called multi-domain. We introduced these two situations because of the slightly differ- 
ent numerical treatment of the problem. Indeed, vacuum regions are treated analytically whereas 
regions filled with plasma are treated numerically. 
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3.1. Single domain 

Let's start with a column density in contact with the inner and the outer conducting wall. The 
motion of the column of plasma is governed by the conservation of charge, the Maxwell-Poisson 
equation and the electric drift approximation, respectively : 



§+div(pv) = 
ot 

+ — =0 

So 

E aB 



(15) 
(16) 

(17) 

E = -V0 (18) 

We recall that the magnetic field remains constant in time in the diocotron regime. We apply the 
standard linear perturbation theory. Introducing perturbations of physical quantities X like electric 
potential, density and velocity components, by the expansion 



(19) 



the eigenvalue problem for the perturbed electric potential <p(r) is expressed as 



1 _5 
r dr 



dip 



1 d 



(oj - mO.) £() r dr\B 



(20) 



For the purpose of our numerical algorithm, it is more convenient to rewrite it as a generalised 
linear eigenvalue problem as follows 



a>-C,„(<f)) - mQ. £,„{(()) + 

The Laplacian operator Xm for each azimuthal mode m is given by 



j:mi<t>) = 



\_d_ 

r dr 



dr 



and the function q,„ is 



m d 
sor dr\ B 



For the electrospheric plasma, we can eliminate p by Eq. (fT4l i and write 



m d 
r dr 



I d I \ rQ. dB^ 
-—(r^Q.) + 



(21) 



(22) 



(23) 



(24) 



r dr ^ ' B^ dr 

Note that this function does not depend on the intensity of B^ but only on its radial profile. However, 
it depends linearly on the amplitude of the differential rotation, meaning that strong departure from 
corotation leads to strong instability as will be shown in the next section. The eigenvalue problem 
Eq. ( |2TI ) is supplemented by homogeneous boundary conditions such that the perturbed potential 
vanishes at the boundaries, namely 



(l){Wx) = 0(VK2) = 



(25) 



3.2. Muiti-domain decomposition 

When vacuum gaps exist between the walls and the plasma, the discontinuity in the density profile 
introduces a Gibbs phenomenon and therefore drastically decreases the efficiency of our pseudo- 
spectral algorithm. To overcome this difficulty, we decompose the space between the two walls into 
three distinct regions, namely 
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- region I: vacuum space between inner wall and inner boundary of the plasma column, the 
solution for the electric potential is denoted by 0i, defined for W] < r < Ri ; 

- region II: the plasma column itself located between Ri and R2, solution denoted by (pn, defined 
forT^i <r<R2 \ 

- region III: vacuum space between outer boundary of the plasma column and the outer wall, 
solution denoted by (pm, defined for R2 < r < W2. 

In region I and III, the vacuum solutions which satisfy the required boundary conditions {4>i{W\) - 
and <^ui(W2) - 0) are given by the usual solutions to Laplace equation in the form 



hir) = A 



,(r) = B 



y.m+ 1 



(26) 



(27) 



where A and B are two constants to be determined by matching the boundary conditions at the 
interfaces between plasma and vacuum. Because no surface charges accumulate on these interfaces 
(the equilibrium density profile should vanish at R\ and R2), the electric field is continuous in the 
whole space. The matching conditions are therefore, continuity of (p and of its first derivatives, 
namely 



*ii(-Ri) 
b[i(Ri) 
MRi) 
b[i(R2) 



(28) 
(29) 
(30) 
(31) 



We eliminate A and B from these equations to get the boundary conditions for 0n- We find that in 
the plasma column, the potential has to satisfy mixed boundary conditions of Robin type 



a\/2 (/>u(Rl/2) + P\I2 0u(^l/2) = 

with the coefficients given by 



a\l2 - 



niR'l'n^ +(m+l) 



P112 = - 



M/2 



R'?n - 



^^\I2 
R\I2 



W2m+1 \ 
^^\I2 



M/2 



om+l 
^^1/2 J 



(32) 

(33) 
(34) 



The boundary conditions Eq. (|32| | are the generalisation of the homogeneous Dirichlet case Eq. 
when vacuum gaps are present. If the vacuum regions were removed such that W\ - R\ and W2 - 
R2, the coefficients /3i/2 in Eq. (|32] | would vanish whereas ai/2 + 0. We therefore find again the 
homogeneous boundary conditions of the previous subsection, Eq. i 



4. NUMERICAL PROCEDURE 

4.7. Single domain 

In the single domain decomposition, the plasma is in contact with both, the inner and outer wall. 
Thus Wi = Ri and W2 - R2- There is no need to distinguish between R and W. The generalised lin- 
ear eigenvalue problem Eq. ( 1211 1 is solved numerically by using a pseudo-spectral Galerkin method 
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taking advantage of the homogeneous boundary conditions imposed on the perturbed electric po- 
tential, (piR\/2) = 0. First, the independent variable r in Eq. (ISTT l is transformed into a new inde- 
pendent variable x such that r € [R\,R2] is mapped into the interval x e [-1,1]. This coordinate 
transformation reads : 
R2-R1 R2+R1 

r^^x^^- (35) 

Second, the unknown eigenfunction is expanded into a set of basis functions {il/k)k>2 (Bovd l200TT l 
defined for « > 1 by : 

4l2n{x) = T2„(x) - 1 (36) 
Ip2n+l(x) = T2„+l(x)-X (37) 

Therefore the boundary conditions are automatically satisfied for each function of this particular 
basis, il)k(+l)\k>2 - 0. Here Tii{x) = cos(A: arccos(x)) are the Tchebyshev polynomials. Let being 
the number of collocation points (jC/Oos/sw i for the Tchebyshev expansion : 



Xi = cos 



.■ = 0,1,..,A^-1 (38) 
The unknown function is expanded according to : 

N-l 

cf>{x)^Yj<l>tMx) (39) 

k=2 

The matrix discretisation of Eq. (|2TJ leads to the generalised linear eigenvalue problem 

L4> = a)M4> (40) 

where we have introduced the unknown vector ^ = {<p2, <pi, (pN-i) and the square matrices are 
defined by : 

Li-i,k-2 = mQ(xi)£„(i//k)(xi) + qm(xi)i//k(xi) (41) 

M,_i,^2 = -Cm(>/^k)(Xi) (42) 

i = l,2,-,N-2 
k = 2,3,..., A?- 1 

The system Eq. ( |40l i is then solved by standard routines to compute the eigenvalues and the eigen- 
vectors. 

4.2. Multi-domain decomposition 

In the multi-domain decomposition, the Galerkin basis Eq. (|36]|-(|37]| is not adapted anymore. The 
new basis functions have to satisfy the following Robin boundary conditions : 

ai.At(-l)+j6i'Al.(-l) = (43) 
a2W+l)+y62-Al.(+l) = (44) 

We look for basis functions expressed by a three term expansion in Tchebyshev polynomials, for 
k>2, like 



tffk(x) = Tk(x) + at Tk-i(x) + bk Tk-2(x) 



(45) 
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For ^k{x) to satisfy the conditions Eq. d43]l-(l44li. the unknowns and bk have to be solutions of 
the following system 

ak[l5x{k-\f-ai]+bk[ai-lix{k-2f] = -ai+Pik^ (46) 

ak [Pi ik-lf + 02] + bk [ai +P2{k- if] = -^2 - P2 k^ (47) 

The procedure to discretise the eigenvalue problem is the same as the one described in the previous 
section. We need only to replace the old basis Eq. ([36]l-(l37Ti by the new basis Eq. ( |45] l. 



5. RESULTS 

Using the aforementioned numerical algorithm, we compute the eigenfunctions and eigenspectra 
of the diocotron instability for various equilibrium density profiles and charges on the inner wall 
for laboratory plasmas. Sec. 15. H and several rotation profiles for pulsar's electrosphere. Sec. 15.21 

In all computations, the external magnetic field follows a decreasing power law with index a > 
such that 



5. 1 . Plasma column 

First, we consider the laboratory plasma confined by some external experimental electromagnetic 
devices. The magnetic field and the density profile are specified as initial data. The rotation profile 
is then deduced from the electric drift approximation Eq. ( fTTT i. 



5.1 .1 . Constant density profile 

The simplest case corresponds to a uniform magnetic field Bq, a - and a constant charge den- 
sity po extending from the origin r = to the outer wall, so Wi - R\ = 0, no inner conducting wall, 
Wi - Ri and also Q-Q (single domain). In this particular case, the equilibrium speed is constant : 

" = (49) 

The functions are all identically zero and the analytical solution to the eigenvalue problem 
Eq. ( |40] | is simply : 

w = m Q (50) 



There is no instability. Our numerical algorithm is able to reproduce the eigenvalues Eq. ( 1501 ) with 
great accuracy taking only a few discretisation points (9 points are already enough). The conver- 
gence is very fast due to the evanescent truncation error introduced by the Tchebyshev expansion. 
The precision achieves easily 10 digits. 



5.1.2. Decreasing density profile 

Next we consider a monotonely decreasing density profile starting from si r - W\ - R\ and 
vanishing ai r - W2 - R2 (single domain) such that 



10 
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The corresponding rotation profile can be found analytically. The electric field induced by the 
plasma is deduced from Eq. ^ whereas the contribution from the wall is given by Eq. (|2]i. Using 
the drift approximation Eq. ( fTTT) . the total radial electric field expressed in Eq. Q and the power 
law index for the magnetic field Eq. ( |48] ), the rotation speed of the plasma column in the most 
general case (a and Q arbitrary) is equal to 



n = — 



eoBo(Rj-Ri) 



1 



Q 



IjTEor^-" Bq 



(52) 



When the plasma is in contact with the inner as well as with the outer conducting wall and the 
applied magnetic field is uniform, there is no diocotron instability, Briggs et al. d 19701 1. Indeed, for 
a = 0, all the eigenvalues we found are real as expected. For the first azimuthal modes m < 19, 
the map of the eigenvalues in the complex plane is shown in Fig.[T]for an uniform magnetic field, 
a = 0, zero charge, Q = Q, Wi = R] = I and W2 = R2 = 20. We normalised the frequencies to the 
rotation frequency co^ = a>^/4a>c - 1/400. Therefore, for each mode, the spectrum is continuous 
and delimited by two integer values, namely 



> a),nl<jJ-[ > -2 m 



(53) 



A similar conclusion applies also when the inner wall is removed (Wi - R\ - 0), see Fig.|2l In this 
case the continuous spectra satisfy 



■ m > a>„Ja>r > —2 m 



(54) 




1 -J5 
a 

S -20 



-23 r 
-30 

-35 ; , ; 

ti 3.3 9 11} 12,3 IS J7.t 

Fig. 1. Eigenvalues for the monotonely decreasing density profile Eq. ( fsTl ). The boundaries are set 
to Wi - R\ - 1 and W2 - R2 - 20. The other parameters are a - and Q - Q. The continuous 
spectrum is shown for the azimuthal modes m - 0..19. The eigenfrequencies have been normalised 
to the rotation frequency cj^ - co^jAoJc - 1/400. 



Finally when switching from an uniform magnetic field to a decreasing power law Eq. ( |48T l with 
for instance a - \, the results are qualitatively unchanged, see Fig. [3] As can be seen by comparing 
the two maps on Fig. [T]and Fig. [3] the radial profile of the magnetic field does not influence much 
the eigenvalues. There is mainly a difiference in scale because the rotation profile for a - I leads to 
higher rotation rate due to the extra factor in Eq. ( |52] l. 
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Fig. 2. Same as Fig.[T]but the inner wall has been removed, i.e. Wi - Ri = 0. 

BCgEn-.'aI.ues 




-doa 



Fig. 3. Same as Fig.[T]but with a decreasing external magnetic field, a - I. Note the difference in 
scale for both cases. 



5.1.3. Gaussian density profile 

Next we consider a situation in which the plasma is not in contact with neither the inner nor the 
outer wall and so Wi < Ri and W2 > R2, (multi-domain decomposition). We choose a Gaussian 
density profile given by 



P^Poe 



-{r-afila-^ 



(55) 



cr is a measure of the profile radial spread while a specifies the location of the centre of the plasma 
column. The associated rotation rate of the plasma column in the most general case {a and Q 
arbitrary) is found analytically to be : 



Q 



Ri-a] 



e ' ' — e 



V2 



erf 



V2 



Q 



(56) 



7.n eqv^ "Bo 

er{(x) is the error function, Abramowitz and Stegun (1 19651 1. In this case, we expect a strong dio- 
cotron instability with growth rates of the same order of magnitude as the rotation speed of the 
perturbation. 

Indeed, taking boundaries such that Wi = I, W2 = 20, Ri - 5 and /?2 = 10 and an uniform 
magnetic field a = 0, no charge on the inner wall Q - 0, the map of the eigenvalues in the 
complex plane is given by Fig. |4] Only the low azimuthal unstable modes are excited, from m = 1 
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to m = 7, the fastest mode being m - A with y,mix ~ 1-15 x 10"^ which is of the same order of 
magnitude as the real part, Re{a>)/m. Removing the inner wall has little effect on the diocotron 




Fig. 4. Eigenvalues for the Gaussian density profile Eq. dSSl l with cr - I and a -1 . The boundaries 
are set to Wi = 1, W2 = 20, Ri =5,R2^ 10 and a = 0. 

spectrum, compare Fig.|4]and Fig.|5] Indeed, the instability is initiated by the presence of a hole in 
the middle of the plasma column. Placing or removing an inner wall cannot suppress or generate 
the instability in this configuration. It can at most affect the rate at which the instability grows. Also, 



C-ldp of ei.£^en'.'.fiiue> = 




Fig. 5. Same as Fig.|4]but the inner wall has been removed, i.e. Wi - 0. 



a decreasing magnetic field a - I does not affect the results in a dramatic way, compare Fig.|4]and 
Fig. |6] The number of excited modes is directly related to the radial spread of the plasma column. 
More precisely, the thiner the plasma ring, the larger the number of unstable modes. In other words, 
decreasing cr will generate more unstable modes. In Fig.|7] we compare three Gaussian distributions 
with different spreads, namely cr = 1, 0.6, 0.4. The cr = 1 profile has 7 unstable modes whereas 
the cr = 0.6 profile has 10 unstable modes and the cr = 0.4 profile 14. Meanwhile, the highest 
growth rate is also increasing significantly. 

Finally, adding a charge to the inner wall does not affect drastically the map of eigenfrequencies, 
see Fig. |8] More unstable modes appear, nevertheless the maximum growth rates remain approx- 
imately in the same order of magnitude. Now that the algorithm has been check and employed 
for different lab plasmas, we study in the next subsection the interesting case of an electrospheric 
plasma. 
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Fig. 6. Same as Fig.|4]but with a decreasing external magnetic field a - I. 
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Fig. 7. Comparison of the highest growth rates for different spreads of the annular plasma column, 
cr = 1 in red 'stars', cr - 0.6 in green 'x' and cr - 0.4 in blue '+'. 
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Fig. 8. Same as Fig.|4]but with a charge per unit length Q - 10. 



5.2. Electrosphere 

The electrospheric non-neutral plasma, as aheady proved in works by Rrause-Polstorff and 
Michel (|1985al l and Petri et al. ( |2002al l. is confined by the rotating magnetised neutron star The 
most important feature is not the density profile but the rotation speed, although both are related in 
an unique manner, Eq. (fT4l i. Plasma shearing motion between different magnetic surfaces leads to a 
kind of Kelvin-Helmholtz instability, at least in the linear regime of the instability (the eigenvalue 
problems look very similar). (We always keep a single domain decomposition such that Wi = Ri 
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and - R-i). Therefore, for electrospheric plasmas, we choose a rotation profile in the plasma 
column that mimics the rotation curve obtained in the 3D electrosphere. To study the influence of 
the profile, we took three different analytical expressions for the radial dependence of Q by mainly 
varying the gradient in diff'erential shear as follows 

Q(r) = 2 + tanh[a (r - ro)] e"^'' (57) 
The values used are listed in Table ([T]l. 



Q. 


a 






3.0 


5 X 10-5 6.0 




1.0 


5 X 10^5 6.0 


nj 


0.3 


5 X 10-5 10.0 



Table 1. Parameters for the three rotation profiles used to mimics the azimuthal velocity of the 
plasma in the electrospheric disk. 

The angular velocity starts from corotation with the star Q = = 1 (normalised to unity for 
convenience) followed by a sharp increase around r = 6 for Qi,2 and a less pronounced gradient 
around r = 10 for Q3. Finally the rotation rate asymptotes twice the neutron star rotation speed, 
Fig.|9] 



3.79 
l.S 

% ^ 

l.i 

1 . 

Z,S t 7.9 10 12.5 iS 17. a 

Fig. 9. Three choices of differential rotation curves in the plasma column for the cylindrical pulsar 
electrosphere, Qi in red, Q.2 in green and Q.3 in blue. 

5.2.1. Uniform magnetic field 

First consider an uniform magnetic field, a = 0. The maximum growth rates for the three rotation 
curves for each mode m are shown in Fig.fTOl 

The profile having the steepest gradient possesses the largest number of excited unstable modes 
because it corresponds to the case where the smallest scales appear, i.e. m large, red stars in Fig.fTOl 
The largest growth rate, for m = 8 its value is ymax = 3.2, is even larger than the maximum rotation 
speed of the plasma column about Umax = 2.9. The diocotron instability operates on a very short 
timescale, comparable to the period of the pulsar. 
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The second steepest profile possesses less unstable modes as we would expect due to the fact 
that only larger scale structures can emerge with this slope of the differential rotation, green 'x' 
in Fig. [TOl The third smooth profile has only four unstable modes, blue '+' in Fig. [TO] We give 




Fig. 10. Fastest growth rate jmax for each azimuthal mode m and the three rotation profile shown 
in Fig.|9]in an uniform external magnetic field, a = 0, Qi in red 'stars', Q.2 in green 'x' and Q3 in 
blue'+'. 



an example of eigenfunctions for the perturbed electric potential in Fig. [TT] It corresponds to the 
fastest unstable mode for the profile Q2. The boundary conditions impose (p{R\/2) = as seen on 
the plot. The corotation radius defined by 
Re(LLi) 



Q(re) = 

is also shown on this plot, depicted by a vertical bar. 

^ igeji fun c't ±.<An = 

3 

I 

^ 

-1 

-3 



(58) 
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Fig. 11. Real (red curve) and imaginary (blue curve) part of the fastest growing eigenfunction for 
Q.2 in an uniform magnetic field. The vertical bar shows the location of the corotation radius. 



5.2.2. Dipolar magnetic field 

Next consider a dipolar magnetic field, or = 3. The maximum growth rates for the three rotation 
curves for each mode m are shown in Fig. [13] 

We can draw the same conclusions as in the previous case. The largest number of excited unsta- 
ble modes is observed for the steepest profile Q3. Comparing Fig.[T2]with Fig. [10] the value of the 
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maximum growth rate for each mode is not strongly affected by the geometry of the magnetic field. 
Here also, we give an example of eigenfunctions for the perturbed electric potential in Fig. [13] It 




Fig. 12. Fastest growth rate y,,,,, v for each azimuthal mode m and the three rotation profile shown in 
Fig.|9]in a dipolar external magnetic field, a = 3, Qi in red 'stars', ^2 in green 'x' and Q3 in blue 
'+'. 

corresponds again to the fastest unstable mode for the profile ^2- The boundary conditions impose 
0(^1/2) = as seen on the plot. Comparing Fig.[T3]with Fig.[TT] close to the inner boundary, the 
real part of the eigenfunctions look very similar Nevertheless when approaching the outer bound- 
ary, the difference becomes appreciable because of the strong departure from an uniform magnetic 
field. 
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Fig. 13. Real (red curve) and imaginary (blue curve) part of the fastest growing eigenfunction for 
Q.2 in an dipolar magnetic field. The vertical bar shows the location of the corotation radius. 



6. CONCLUSION 

We developed a pseudo-spectral Galerkin code to compute the eigenspectra and eigenfunctions of 
the diocotron instability for arbitrary magnetic field configurations and density profiles. The code 
is very efficient in computing the fastest growing modes. When increasing the number of points, 
convergence is reached quickly, in most cases, less than 200 points are needed. Application to the 
pulsar shows that the diocotron regime gives rise to instabilities with growth rate comparable to 
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the rotation period of the neutron star This confirms the results of our previous work (Petri et 
al l2002bl) and shows that the precise geometry, let it be a flat disk or a column of plasma, does not 
affect drastically the diocotron instability. The required ingredients are only non-neutral plasmas 
and shear velocities. 

In a forthcoming paper, we plan to investigate the full non-linear development of the diocotron 
instability by means of 2D electrostatic PIC simulations. We will also add an external source 
of charge feeding the system and demonstrate that particle transport across the magnetic field 
line is possible. Note that this has already been observed in some experiments by Pasquini and 
Fajans ( |2002] ). 

Extension to plasma moving with relativistic speeds is also envisageable, in particular the gen- 
eralisation to the electromagnetic non-neutral instability regime, the so-called magnetron instabil- 
ity. Last but not least, the influence of finite temperature in the plasma on the diocotron instability 
would require a kinetic treatment of the stability via the Vlasov-Maxwell equation. 
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